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DYNAMICAL SYSTEM ANALYSIS OF REYNOLDS STRESS CLOSURE EQUATIONS 


SHARATH S. GIRIMAJI* 


Abstract. In this paper, we establish the causality between the model coefficients in the standard 
pressure-strain correlation model and the predicted equilibrium states for homogeneous turbulence. We 
accomplish this by performing a comprehensive fixed point analysis of the modeled Reynolds stress and 
dissipation rate equations. The results from this analysis will be very useful for developing improved pressure- 
strain correlation models to yield observed equilibrium behavior. 

Key words, turbulence modeling, pressure-strain modeling 

Subject classification. Fluid Mechanics 

!• Introduction. The equilibrium states of benchmark turbulent homogeneous flows have long been 
used to develop, calibrate and validate pressure- strain correlation models, Speziale et al (1991). The bench- 
mark flows have typically been plane shear and strain-rate dominated flows such as plane strain and ax- 
isymmetric expansion/contraction. When these pressure-strain models are used to compute slightly more 
complex flows which contain the effects of both rotation and strain (elliptic flows), the model results are 
inconsistent with linear stability theory and direct numerical simulation (DNS) data, (Speziale et al 1996, 
Blaisdell and Shariff, 1996). Improved predictive capability of elliptic flows requires a better understanding 
of the physics of turbulence as well as the dynamics of the model equations. This calls for an intimate un- 
derstanding of the model equations, especially the causality between the model constants and the predicted 
long-time behavior. For the special case of shear flow in rotating coordinate frame, the equilibrium state 
predicted by two-equation turbulence models and linear pressure-strain models have been investigated by 
Speziale and Mac Giolla Mhuiris (1989, 1990). These studies have lead to models with better predictive 
capability in those selected flows. 

The objective of the present study is to perform dynamical system analysis, i.e., fixed point and bifurca- 
tion analyses, of the anisotropy, kinetic energy and dissipation closure equations with aid of representation 
theory for all elliptic flows (with two-dimensional mean velocity field) to establish for the first time exact 
(analytical) relationship between turbulence model (pressure-strain correlation and dissipation equation) 
coefficients and the asymptotic behavior of the equation. The results from our study enables us to classify 
completely the turbulent asymptotic state predicted by a model as a function of the mean strain and rotation 
rates. The knowledge of this causality can be used in closure model development to yield the required asymp- 
totic (equilibrium) behavior. Quasilinear pressure- strain model which includes the entire class of Launder, 
Reece and Rodi (1975), referred to as LRR, models are considered here. This study also sheds light on the 
strategies that can be employed to develop non-equilibrium algebraic Reynolds stress models starting form 
the Reynolds stress closure equations. 

In this study we will restrict ourselves to turbulence in inertial frames which includes elliptic flows. 
Turbulence modeling in non- inertial rotating frames will be considered in future work. 

* ICASE, M/S 403, NASA Langley Research Center, Hampton, VA 23681-0001 (Email: girimajWicase.edu ). This research 
was supported by the National Aeronautics and Space Administration under NASA Contract Nos. NAS1-97046 and NASl- 
19480 while the author was in residence at the Institute for Computer Applications in Science and Engineering(ICASE), M/S 
403, NASA Langley Research Center, Hampton, VA, 23681-0001 
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2. Turbulence Closure Equations. In homogeneous turbulence, the exact Reynolds stress transport 
equation in an arbitrary inertial reference frame is given by 


(i) 


du{Uj 


dt* 


— 4*ij • 


The terms, respectively, are the time rate of change, production {Pij), dissipation (e^) and pressure-strain 
correlation (0»j) of Reynolds stress: 


( 2 ) 


p . _ dUj dUi 

*ij — ^ UjUy c , 

ox ) fc dx k 


dui duj dui duj . 

~ dx k dx k ' <l>ij ~ pi dx j + dx?- 


The production and dissipation rate of turbulent kinetic energy are, respectively, P = \Pa and e = ^£u- In 
high Reynolds number flows, dissipation is generally treated as being isotropic: 

(3) e 2 


i £dij - 


Closure models are needed for the pressure-strain correlation (<j>ij) and dissipation rate (e). 

In this study, we focus on the quasilinear class of pressure-strain correlations models of the general form: 


( 4 ) 


<t>ij = -(C?£ + C}P)bij + C 2 KS ij + 

bmnS mn 5ij) + 


c 3 K(b ik s; k + b jk s; k 3 
c 4 K(b ik w; k + b jk w; k ), 

where the C's are numerical constants and 

( 5 ) 


5 . = h™i 

2 { dxj 

_ UiUj _ 

v 2 K 


dUj 


dxi 

W 


)i 


w* 

v 


1 

2 dxj dxi 


We choose this form of the pressure-strain model for two reasons. First, this form of the model permits 
analytical treatment of the asymptotic behavior. Second, this is the form most frequently used in practical 
Reynolds stress closure calculations: this form includes all linear-pressure strain models (e.g. LRR model) 
and some of the non-linear models (such as the quasi- linear SSG model) can also be reduced to this form 
near equilibrium. For the LRR model, the coefficients are 

C? = 3.0; C\ = 0.; C 2 = 0.8; C 3 = 1.75; C 4 = 1.31 
The coefficients for the Gibson and Launder (1978) model are 

C? = 3.6; C} = 0.; C 2 = 0.8; C 3 = 1.2; C 4 = 1.2 
For the quasilinearized SSG model, the coefficients are 

C? = 3.4; C\ = 1.8; C 2 = 0.36; C 3 = 1.25; C 4 = 0.4 


The anisotropy evolution equation can be derived from the Reynolds stress equation and the pressure 
strain correlation model: 

(6) ^ = -bijiL 0 ^ - L\b mn S* mn ) + L ^ 

+L 3 (bikSjk 4- bjkS ik — ~bimSi m Sij) 

+L 4 (b ik w; k + b jk w; k ) 
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The pressure-strain correlation model coefficients are redefined as: 


( 7 ) 


L® = C® — 2; L\ = 2 C{ + 4 ; £2 = C2 — 
7 / 3 = C3 — 2; Z/4 = C4 — 2. 


The turbulent kinetic energy evolves according to 


( 8 ) 



e, 


and the modeled evolution equation of dissipation is 


( 9 ) 


de 

dt* 




The model constants C e 1 and C e 2 are typically given values of 1.44 and 1.90 respectively. 

Equations (6), (8) and ( 9 ) constitute the second order closure equations in homogeneous turbulence. 
These equations can be non-dimensionalized using the norm of the deformation rate tensor: 


( 10 ) 

The non-dimensional quantities are 
( 11 ) 


n = s: j s: J + w* j w* J , 

Sij = w {j = wz/y/tj; 

dt = y/rjdt*; lj = e/( y /rjK) i 


where uj is the relative strain rate, i.e., the ratio of the turbulence to mean flow strain rates. In the 
above equations, asterisk is used to represent dimensional quantities and the corresponding non-dimensional 
quantity is written without the asterisk. In dimensionless time, the anisotropy transport equation is 

( 12 ) = -M£?« - L\b mn S mn ) + L 2 S l3 

2 

-\r L^Sikbkj " 1 “ bikSkj 2 bmn Smn ) 

+Lt(W ik b kj -b ik W kj ). 


Two points are worthy of note here. In dimensionless time the anisotropy evolution is (i) independent of the 
magnitude of deformation (77) and depends only on the ratio of strain to rotation rate; and (ii) dependent only 
on the relative strain rate and not individually on kinetic energy and dissipation. The evolution equation of 
the relative strain rate, w, is easily obtained from those of the turbulent kinetic energy and dissipation: 

(1 3 ) ~ = -MCel - l ) b rnnSmn ~ (C e2 - l)uA 

The first term on the right hand side of equation ( 13 ) represents the production of the relative strain rate 
whereas the second terms represents its destruction. 

In homogeneous turbulence, the above non-linear dynamical system of equations represents an initial 
value problem. The Reynolds stress, kinetic energy and dissipation may grow unbounded from their initial 
values. However, it is known that for some benchmark flows such as (homogeneous) plane shear, plane strain 
and axisymmetric expansion/contraction the normalized turbulence parameters bij and u ; evolve from their 
specified initial conditions according to the equations and asymptote to finite- valued fixed points (equilibrium 
turbulence), provided such a state exists. The fixed point or the equilibrium state of turbulence is described 
by 

/, dba duj 

( 14 ) = 0; and — = 0. 

dt dt 
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2.1. Representation theory. We are interested in the asymptotic behavior, long after the influence 
of the initial conditions has diminished. At this stage of flow evolution, it can be argued that the Reynolds 
stress can be a tensorial function of only the mean strain and rotation rates. Representation theory can 
then be invoked to determine the most general tensor function that can be constructed with the strain 
and rotation rates. The most general, physically permissible tensor representation for the Reynolds stress 
anisotropy in terms of the strain and rotation rates in the case of two-dimensional mean flow is given by 
(Girimaji 1996a) 

(15) 
where, 

(16) r)i = SijSij\ and 772 = 

In the above equations G\ - G 3 are yet to be determined scalar functions of the invariants of strain and 
rotation rate tensors. During the evolution of b t j, G\ G 3 will be functions of time as well. Also note that, 
by definition, 771 + 772 = 1 . 

The representation for b X j is now substituted into equation ( 12 ) and the resulting equation is simplified 
using the following identities valid for all two-dimensional mean flows: 

(17) S ik S kj = ^rnSlf;W ik W kj = 

SikSkiSij — -r)iSij;SikWkiSij = 

W^ik^klWlj = -IftSijiSrnnbrnn = G\7] i, 

where and Sij are the two and three dimensional delta functions respectively. We retain terms up to 
quadratic power (in strain and rotation rate) in their original form and invoke their two-dimensional property 
only when these terms appear in cubic and higher power terms. For example, S xk S k j is retained as such 
when it appears by itself: whereas, when it appears as a part of a cubic or higher power term we invoke 
S ik Skj = to write S xk S k jSji — ~ r}\Su . By invoking the two-dimensional property for reducing only 

cubic and higher power terms, it is hoped that the three-dimensional effect is approximately accounted for 
up to the quadratic term. Using these rules, we write 

(18) Sikbfcj ~h bikSkj — ^bfnfiSfnfiSij = —TJiG^Sij 

+2G\(SikSkj — 

w xk b k3 - b ik W kj = -GxiSikWkj - W ik s kj ) 

-\~2rj2G2Sij . 

Substitution of equation (15) into equation ( 12 ) yields the anisotropy evolution equation in terms of coef- 
ficients G 1 - G 3 . After substitution, the coefficient of each tensor on either side of equation (12) has to be 
equal due to the linear independence of the generators of the integrity basis. Comparing the coefficients of 
the three representation tensors on either side we get 

( 19 ) + G\(L^uj - L\Gir]i) — L 2 + -LzG$r]i 


bij = G ! + G 2 (S ik W kj - W ik S kj ) 
+G 3{Si k S k j — -rji Sij), 
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-f-2Z/4?72^J2 

+ G2(L®C0 — L\G\T]\) = — GiL4, 

+ G^L^uj — L\G\T}i) = 2G r iL3. 

Equations (19) along with (13) constitute the new non-linear system of evolution equations for Reynolds 
stresses in homogeneous turbulence. 


3. Dynamical System Analysis. Fixed point and bifurcation analyses of the new system is now 
performed. 

The fixed point equations (14) can now be restated as 


(20) _ dG2 _ dG$ _duj_ 

dt dt dt dt 

Using the notation that the fixed point values are denoted by a superscript 0, the algebraic fixed point 
relations are (using b mn S mn = G^i): 


(21) 0 = 2 u°(C el - 1)0?*?! + (C e2 _! - l)(u,°) 2 

0 = -C°(L°u, - + L 2 + ^L 3 G° 3m 

+2L4772G2 

0 = Gg(Z,?w - L\G\rn) + G\L a 

0 = G^lo - L\G%) - 2G°L 3 




\uj 


*0 ujulviaa ui ii/auo tu Lilt live liACl 

\ = Oi ^2 + \ L $ Gzti \ + 2L4772G2 = 0 ]; 

/^0 1 r\ /^0 Z/4 /"'rO 2^1 

— n-r,, 1 ^3 — “ 


This system of equations 
[o;° = 0 

0, G i = — 7= . . . 

°’ G ? = +*9i. G ° = ^, 

’° = -2g^G?r, l5 G° = -$L, Gj = G° 

l w ° = G °i = +S> ^ 


£4 /n<) _ 2L 3 

£*171 1 3 L*r/i J 

L4, /^(O 2L 3 I 

°3 - 


C e2 -l^ l7/1 ’ ^ v/m ’ ^2 “ L*rji ’ ^3 “ 

In the above equations Qi and Q* are defined as 


JLLsl] 


(22) Qi = 

Q* = 

where 

(23) 


LJ 


2(^1 


-m 

m 


L2 2 L 2 LI 1-ift 

L* 3 L* 2 L* 2 jji 


L*=2L?§i— 1 + L}. 


■ / c2 


It should be pointed out that a negative value of (7* is consistent with a gradient-diffusion type effect and 
energy flow from the mean to the fluctuating velocity field. A positive value would imply count-gradient 
diffusion and a negative value of production - i.e, flow of energy from the turbulent fluctuations to the mean 
flow. 
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3.1. Bifurcation analysis. The qualitative behavior of the solution of a set of differential equations 
may depend on the parameters of the system. For example, the nature and even the number of fixed points 
of a system can change with changing parameter values (bifurcation). The parameters of the present system 
of equations are the constants of the pressure- strain correlation model and the mean strain-rate to total 
deformation ratio r)\. (Note that, by definition, 772 = 1 — 771, and hence 772 is not considered an independent 
parameter.) In this paper, we restrict our attention to bifurcation due to 771 alone for a given pressure-strain 
correlation model. 

The system has five fixed points when 771 is such that both Q\ and Q* are real. If only one of Q\ or Q* 
is real, then the system has three fixed points. If both are not real, then there is only one fixed point. It is 
important to establish the conditions under which Q\ and Q* are real. 

For Q 1 to be real we require 


( 24 ) 


i ^ ( ^3 \2 

' L\ + 3 ( L\ } 


2( Zr )2 V^- 0, 

^1 vi 


implying that 771 should be greater than a critical value 77“: 


( 25 ) 


m>Vi = 


2 L\ 

-L 2 L\ + lLl + 2 L\' 


Similarly, Q* is real only when 


( 26 ) 


Vi>Vi = 


2L4 

-L 2 L*+ |L| + 2Lf 


For all the pressure-strain correlation models considered in this paper, L 2 is negative and L* > L J leading 
to 


( 27 ) 1 > 77“ > ttf > 0 . 

When 771 is in the interval (1, 77“ ), the system has five fixed points. The system has three fixed points in the 
interval (77® , 77J). Finally, for 771 < 77^, the system has only one fixed point. 

In summary, the nonlinear set of equations governing the (modeled) evolution of the Reynolds stress has 
two bifurcation points 77“ and r}\. The behavior of the solution for various values of 771 depends upon the 
stability of the various fixed points. 

3.2. Stability of fixed points. In order to establish the stability of a fixed point, we need to examine 
if any small perturbation of the system away from the fixed point eventually returns to the fixed point after 
a sufficiently long time. The most expeditious way of establishing this is by investigating the Jacobian of 
the system at the fixed point. If an eigenvalue of the Jacobian is negative, solution trajectories are attracted 
towards the fixed point along the corresponding eigenvector. On the other hand, if an eigenvalue is positive, 
the solution trajectory is repelled away from the fixed point along the corresponding eigenvector. For a fixed 
point to be stable all the eigenvalues must be negative, so that all trajectories in the neighborhood of the 
fixed point are attracted towards the fixed point. If all permissible initial conditions are attracted to the 
fixed point, such a fixed point is called globally asymptotically stable. The set of all initial conditions that 
ultimately evolve to a stable fixed point is called the basin of attraction of that fixed point. For a nonlinear 
set of equations, such as the one considered here, it is difficult to establish the basin of attraction and will 
not be attempted here. We will only seek to establish the local asymptotic stability of each of the fixed 
points. First, the various types of fixed points are fisted. 
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1. When all the eigenvalues of the Jacobian are real and negative, then the fixed point is a stable fixed 
point also called an attractor or a sink. 

2. If all the eigenvalues are real and positive, then the fixed point is a source or repellor. All solution 
trajectories in the neighborhood of the fixed point are repelled away from it. 

3 . If the eigenvalues are real with some positive and the others negative, then the fixed point is of the 
saddle type. The solution trajectories are attracted towards the fixed point in some directions and 
repelled away in other directions. 

4 . If the eigenvalues are complex and the real part is positive, then the fixed point is a spiral source. 
The solution to the system fluctuates about the fixed point with the amplitude of the fluctuation 
getting larger with time. 

5 . If the eigenvalues are complex with a negative real part, then the fixed point is a spiral sink. The 
solution is oscillatory with decreasing magnitude about the fixed point. 

6. If the eigenvalues are purely imaginary, then the fixed point is classified as a center. The asymptotic 
solution then displays an oscillatory behavior. 

The sink and spiral sink fixed points are attracting fixed points and, hence, stable. Saddle and source fixed 
points are unstable. 

Fixed point # 1 .. This fixed point is given by 

u ; 0 = 0, G i =0, L>2 ~ T 3 G 3771 -T 2 L 4772 G 2 = 0. 

This is the only fixed point that exists for the entire range of 771 values. 

The Jacobian at this fixed point is 

0 0 0 0 

0 0 2X4(1-17!) \Lzrn 

0 L\t}\G*2 — ^4 0 0 

0 L\r}iG® + 2 L 3 0 0 

The eigenvalues of the Jacobian are 

*1=0, *2 = 0, *3 = +y/rjiL\Q, *4 = ~y/rjiL\Q, 

and the corresponding eigenvectors are 



v a = [1,0, 0,0], 

v 3 = [ 0 , + y/rj\L\Q , L\t]iG2 — T4, L\rj\G^ + 2L3], 
v 4 = [ 0 , —y/rjiL\Q, L\t)\G 2 — L4, L\t)\G 3 + 2L3]. 


Since two of the eigenvalues are zero, this fixed point is classified as a non-hyperbolic fixed point. So long 
as Q is real (77 1 > 77J), *3 is positive and *4 is negative leading to this fixed point being a saddle and hence 
unstable. When 771 < 77^ the eigenvalues *3 and *4 are both purely imaginary; the fixed point is a center and 
the asymptotic solution is oscillatory. The stability of non-hyperbolic center fixed point cannot be gleaned 
from a linear approximation. The stability can be evaluated only from a center manifold theory which is 
not attempted here. It suffices here to say that for 771 < 77J, the long-time behavior of the Reynolds stress in 
oscillatory. 
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Fixed points # 2 and # 3 .. These fixed points are given by 


u> 


= 0. g» = ±-U>„ <4 = -4t, g$--t=T’ 

v/^1 ^ »7i L Vi 


and exist only in the range 771 > 77 ?. The Jacobian at this fixed point is 

/ 0 0 0 0 \ 

— G° 2G* 2L4(l-m) |L 3 77 i 

— 0 G** 0 

0 G** 


(30) 


V -G? 0 


where G* = G^L} + C e i - 1) and G** = G 0 im (L\ + 2G e i - 2). 
The eigenvalues are given by 


(31) 


Ai = -2(G e i - 1)G°77 1 , A 2 = 2L}G? m , 
A3 = A 4 = L{GiT)\. 


The eigenvectors are, 


(32) 


v 1 = [^(Ll + 2G el -2), 

G ! i r ?i(-^'i + 2 G e2 — 2 ) + 2L 4 »72 + 5^3771 
2G°7 /1 (Ll+G e i -1) 
v 2 = [0, 1,0,0], 

v — (0 2L 4 (1 — 771 ) 

V3 - [0 ’“ G°L\m ’ ’ 


, 1 , 1 ], 


v 4 = [0, 


3G,L] 


, 0 , 1 ]. 


Since all the eigenvalues are non-zero this an hyperbolic fixed point. In all the models considered, L\ and 
(G e i — 1) are both positive and 771 is positive by definition. As a result, irrespective of the sign of G, . some 
of the eigenvalues will be positive and others negative and the fixed point is a saddle. These two fixed points 
are unstable when 771 > 77 “ and do not exist otherwise. They do not play an important role in the long-time 
behavior of the Reynolds stresses. 

Fixed point This fixed point is given by 


(33) 


w° = - 2 ^— g; 

e2 “ 1 

n 0 _ ^4 _ 2 L 3 

^9 — — 1 ^3 — — “ 




L*rh 


L* 771 


This fixed point exists only for 771 > r}\. Due to the complex nature of the Jacobian, it is difficult to obtain 
all the eigenvalues and eigenvectors symbolically. However, one eigenvalue is easily obtained by inspection: 

(34) A, = G° 77iL*; v, = [0, 0, 1, -6 ^—^ ]. 

The eigenvalues evaluated numerically are plotted in Figure 1 (a) as a function of the parameter 771 for 
the linearized-SSG pressure-strain correlation model. All of the eigenvalues are non-zero (hyperbolic fixed 
point) and negative indicating that this is an attractor. (Note that the quantity plotted is the negative of 
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the actual eigenvalues.) Another important point to be gleaned from the the figure is that eigenvalue A4 is 
always about an order of magnitude smaller than the other eigenvalues. This indicates that the evolution 
equations evolve slowly along the eigenvector associated with A 4 and rapidly along all other direction. The 
eigenvector direction corresponding to A 4 is shown in Figure 1(b) and it is almost coincident with the u> axis. 
The behavior of the eigenvalues of this fixed point with LRR pressure-strain correlation model is qualitatively 
and, even, quantitatively similar to that of SSG model. 

Fixed point #5 . . This fixed point 


(35) 


= -2 


C, 


el 


1 




^4 2^3 


Vm 


m 


g° 3 




also exists only for rji > r)\. The eigenvalues calculated numerically are all positive indicating that the fixed 
point is a source or a repellor and, hence, unstable. 

All the stability and bifurcation results for Gh are summarized in Figure 2. The bifurcation diagram of 
to is given in Figure 3. 


4. Summary and Discussion. The long-time behavior of the standard Reynolds stress closure equa- 
tion depends upon the value of the parameter rji , the fraction of deformation that is strain. The asymptotic 
behavior undergoes bifurcation at T ) i = rft. For values of 771 higher than ry, , the solution of the equation set 
asymptotes monotonically to an attractor. For smaller values, the long-time behavior is oscillatory. 

Asymptotic behavior for rji > rf\.. The solution is attracted to fixed point 4. The production to dissipa- 
tion ratio at this equilibrium state is 


(36) 


— (equilibrium) 


Ce2~ 1 
C e i-\ 


independent of r )\ . In appropriately normalized time, the turbulent kinetic energy and dissipation grow 
exponentially at rates independent of rji : 


(37) 


dlnK 

dr 



dine P 

^ r _c el7 -c e2 , 


where r = u>ot. 

Longtime behavior for rji < r}\.. The long-time behavior in this case is dictated by fixed point 1 which 
is non- hyperbolic (has zero eigenvalues). Linear analysis about the fixed point is inadequate to determine 
its stability. Numerical calculations indicate that (i) to decays monotonically (exponentially?) to zero, (ii) 
G\ oscillates about zero, (iii) G 2 and G3 converge to non-zero values, (iv) with decreasing rji the asymptotic 
growth rates of kinetic energy and dissipation decrease. In fact, for small enough rj\ values kinetic energy 
and dissipation decay in time leading to relaminarization. 

Implications to pres sure- strain correlation modeling.. It is reasonable to demand that future models 
yield bifurcation diagrams consistent with experimentally (laboratory and numerical) observed behavior of 
the Navier-Stokes equation. At the very least, for the sake of qualitative consistency with true physics, the 
bifurcation points of the model equation set should coincide with that of Navier-Stokes statistics. 

The DNS data of Blaisdell and Shariff (1996) appears to indicate that the asymptotic growth rate of 
kinetic energy is always positive and nearly independent of rji : 


(38) 


dlnK 

dt 


~ Const. (> 0), for 1 > 771 > 0, 
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a result that is consistent with linear stability analysis. The lack of qualitative (or, even, quantitative) change 
in the behavior of the kinetic energy growth rate appears to imply that there is no bifurcation, contrary 
to model predictions. This can be interpreted as the bifurcation taking place at 77 1 = 0. If this is indeed 
true (further work is currently underway), then, in order for the pressure strain model to be consistent with 
physics we need: 


(39) 




2L\ 


0 , 


—L2L* + + 2 T4 

leading to 1/4 = 0, implying C 4 = 2. 


This should be the behavior of C 4 in the limit of vanishing values of 771 . The currently used values of C 4 
yield good agreement with data in homogeneous shear and strain dominated cases. This clearly suggests 
that C 4 and perhaps all other coefficients should be functions of 771 and 772 . These ideas will be explored in 
future works. 

Non- equilibrium algebraic stress modeling .. Non-equilibrium algebraic Reynolds stress model can be 
considered as the approximate solution of the Reynolds stress closure equations away from the equilibrium 
state. When 771 > 77J, the solution approaches the equilibrium state (fixed point 4) along the invariant 
manifold corresponding to eigenvector V 4 for nearly all initial conditions. This is because the corresponding 
eigenvalue (A 4 ) is much smaller than the rest (Figure la). (From a random initial condition, the solution 
evolution along the other eigenvector directions is very rapid, so that the trajectory after a short initial stage 
is nearly aligned with the slowest eigen-direction.) In Girimaji (1996b), it was suggested that the invariant 
manifold of V 4 be used as the non- equilibrium Reynolds stress model and a strategy for determining the 
manifold was presented. When 771 < 77J , the strategy for non-equilibrium algebraic modeling in not yet clear. 
Central manifold reduction may be the answer and that approach is currently under investigation. 
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Eigenvalues of stable fixed point (SSG) 



Eigenvector Components of X 4 



Normalized Strainrate rj x 

Fig. 1 . (a) Eigenvalues of fixed point 4 as a function of tj . (b) Components of eigenvector corresponding to eigenvalue 

A 4 . 
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Bifurcation diagram of G x (SSG) 



0.0 0.2 0.4 0.6 0.8 1.0 

Vi 

Fig. 2 . Bifurcation diagram of G\. Fixed point <?}: unstable (saddle) for 77 > r)\ and center for 77 < rft. G 2 , G 3 
unstable (saddle). G A stable (attractor). G*>. unstable (repellor). 


Bifurcation diagram of cj (SSG) 



0.0 0.2 0.4 0.6 0.8 1.0 

Vi 

Fig. 3 . Bifurcation diagram of u ;. Fixed point lj 1 (— a; 2 = cu 3 ) is unstable (saddle) for 77 > 7]l\ and stable for 77 < 771. 
o> 4 (= u; 5 ) is stable (attractor). 
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